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Abstract 

This  report  reviews  the  theory  and  practical  estimation  of  condition  numbers  for  the 
nonsymmetric  eigenvalue  problem.  The  report  provides  a  manual  for  using  LAPACK 
subroutines  STRSNA  and  STRSEN  to  estimate  condition  numbers  for  individual  eigenvalues 
and  eigenvectors,  multiple  (or  clustered)  eigenvalues,  and  invariant  subspaces. 
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1      Introduction 

We  review  the  theory  of  condition  numbers  for  the  nonsymmetric  eigenproblem,  and  describe 
algorithms  for  estimating  them.  We  provide  a  manual  for  the  LAPACK  subroutines  STRSNA 
and  STRSEN,  which  compute  these  condition  numbers  for  matrices  in  Schur  canonical  form. 
We  assume  the  reader  is  familiar  with  the  basic  theory  of  the  nonsymmetric  eigenproblem: 
eigenvalues,  right  and  left  eigenvectors,  multiple  eigenvalues  and  right  and  left  invariant 
subspaces. 

The  condition  number  of  a  problem  measures  the  sensitivity  of  the  solution  to  small 
changes  in  the  input.  We  call  the  problem  ill-conditioned  if  its  condition  number  is  large, 
and  ill-posed  if  its  condition  number  is  infinite.  We  may  use  condition  numbers  to  bound 
errors  in  computed  solutions  of  numerical  problems. 

We  illustrate  this  with  a  simple  example.  It  is  well  known  that  the  condition  number  for 
solving  a  system  of  linear  equations  is  k{A)  =  ||/1||-||.4~^||,  where  ||-||  is  any  matrix  operator 
norm  (we  will  be  more  specific  about  norms  later).  Suppose  that  linear  system  Ax  =  b  is 
solved  via  Gaussian  elimination  with  partial  pivoting,  or  some  other  stable  scheme.  Let  x 
be  the  computed  solution.  Then  one  may  bound  the  error  by: 

\\x  -  x\\ 

"  „    „  "  =  0{macheps)  ■  k{A) 

where  macheps  is  the  machine  precision.  The  size  of  the  constant  implicit  in  the  O(-) 
notation  depends  on  the  size  of  the  matrix,  pivot  growth,  etc. 

Condition  numbers  may  be  expensive  to  compute  exactly.  For  example,  computing 
k{A)  for  even  the  simplest  matrix  norms  is  three  times  as  expensive  as  solving  Ax  =  b  in 
the  first  place.  Therefore,  one  usually  uses  an  inexpensive  estimate  in  place  of  the  exact 
k{A).  For  example,  a  method  for  estimating  k{A)  is  included  in  LINPACK,  which  costs 
just  O(n^)  extra  beyond  the  0{n^)  cost  of  solving  Ax  =  b.  The  price  one  pays  for  using  an 
estimate  is  occasional  (but  hopefully  rare)  misestimates  of  k{A).  Years  of  experience  with 
the  LINPACK  estimator  attest  to  its  reliability,  although  examples  do  exist  where  it  can 
underestimate  k{A)  badly. 

The  codes  we  discuss  for  the  nonsymmetric  eigenproblem  will  also  use  such  condition 
estimators.  Here,  the  savings  will  be  even  greater  than  for  linear  equation  solving:  an  0{n'^) 
estimator  using  0{n)  workspace  in  place  of  an  0{n^)  exact  solution  using  O(n^)  workspace 
in  some  cases. 

Our  condition  estimators  wiU  compute  two  quantities,  the  reciprocal  of  a  condition 
number  for  an  eigenvalue  (or  cluster  of  eigenvalues),  and  the  reciprocal  of  a  condition 
number  for  an  eigenvector  (or  invariant  subspace).  We  compute  reciprocals  of  condition 
numbers  to  avoid  overflow;  an  infinite  or  overflowed  condition  number  is  indicated  by  a 
zero  reciprocal.  By  combining  these  two  values  in  simple  algebraic  formulas,  a  great  deal 
of  detailed  information  about  the  eigenproblem  can  be  obtained.  This  report  will  describe 
both  these  basic  condition  numbers  and  these  formulas. 

Our  condition  numbers  will  measure  the  changes  in  the  eigenvalues,  right  eigenvectors, 
means  of  clusters  of  eigenvalues,  and  right  invariant  subspaces  of  a  matrix  A  when  a  pertur- 
bation E  is  added  to  it;  our  bounds  will  be  functions  (usually  multiples)  of  \\E\\.  This  may 
be  used  to  estimate  the  error  in  solutions  computed  by  LAPACK  routines  because  they  are 


backward  stable,  i.e.  they  compute  the  exact  eigendecomposition  of  a  matrix  A  +  E  where 
A  is  the  input  matrix,  and  H^H  =  0{Tnacheps)\\A\\.  We  measure  changes  in  eigen\-ectors 
and  invariant  subspaces  by  their  change  in  angle;  we  discuss  the  angle  between  subspaces 
in  more  detail  in  section  5.  Our  condition  numbers  yield  both  asymptotic  bounds,  which  are 
accurate  only  w^hen  the  norm  ||£|i  is  small,  and  global  bounds,  which  work  for  all  ||£||  up  to 
a  certadn  upper  bound,  whose  size  depends  on  the  problem  and  may  be  large  or  small.  We 
show  how  to  obtain  these  upper  bounds  on  ||£||  as  well. 

We  illustrate  the  reason  for  providing  such  a  variety  of  bounds  with  an  example.    Let 
A,,  be  11  by  11  of  the  following  form: 

0      1  0 


^v  —  ■  • .     1 

0     0 

.5 

Here,  blank  entries  are  also  zero.  Thus,  Arj  is  a  block  diagonal  matrix  with  a  10  by  10  block 
at  the  upper  left  and  a  1  by  1  block  at  the  lower  right.  When  rj  =  0,  the  upper  left  block  is 
a  10  by  10  Jordan  block  with  a  single  multiple  eigenvalue  at  0  and  a  single  right  eigenvector 
V  =  [1,0, ..  .,0]-^.  Such  a  matri.x  is  called  defective.  For  small  nonzero  t]  the  eigenvalues 
become  distinct  numbers  all  with  absolute  value  t]-^,  and  eigenvectors  which  have  rotated 
away  from  v  by  about  77-^  radians.  When  77  =  10~'°.  tj-^  =  .1,  a  much  larger  change.  In 
this  case  we  call  the  eigenvalue  at  0  and  its  associated  eigenvector  ill-posed,  because  their 
sensitivity  is  not  proportional  to  the  norm  of  the  perturbation  77,  but  a  root  of  tj. 

The  practical  solution  to  this  problerv  is  to  consider  this  matrix  as  having  a  cluster 
of  10  eigenvalues  near  zero  with  a  single  invariant  subspace  which  is  spanned  by  all  their 
eigenvectors,  as  well  as  a  single  eigenvalue  near  .5  with  its  eigenvector.  The  mean  of  this 
cluster  of  10  eigenvalues  will  be  much  less  sensitive  to  small  perturbations  than  the  individ- 
ual eigenvalues  (in  fact  it  will  be  independent  of  77  in  this  example).  For  small  enough  ||£^||, 
our  asymptotic  error  bounds  will  show  that  the  mean  of  the  cluster  of  10  eigenvalues  near 
0  of  ylo  +  -E^  is  bounded  by  ||£^||  (see  Bound  4  below);  i.e.  the  mean  is  very  well-conditioned. 
Similarly,  the  eigenvalue  near  .5  can  also  only  change  by  H^H  for  small  enough  ||£||  (Bound 
2).  The  invariant  subspaces  will  also  be  much  less  sensitive  than  the  individual  eigenvectors. 
In  this  example,  the  right  invariant  subspace  belonging  to  the  cluster  of  10  eigenvalues  near 
0  is  spanned  by  the  first  10  columns  of  the  11  by  11  identity  matrix  independent  of  77;  more 
genercdly  our  bounds  will  say  that  for  small  H^H  the  right  invariant  subspace  can  rotate 
by  at  most  2731||£||  radians  (Bound  6).  The  eigenvector  belonging  to  the  eigenvalue  .5  is 
equally  insensitive  in  this  example. 

This  illustrates  our  asymptotic  error  bounds,  valid  for  sufficiently  small  \\E\\.  In  contrast, 
our  global  bounds  give  bounds  valid  for  all  ||£||  up  to  an  upper  bound  which  we  also 
estimate.  The  matrix  An  illustrates  the  source  of  these  upper  bounds  on  ||£'||.  Suppose  we 
make  77  =  2~^°  ^  .001;  then  one  of  the  eigenvalues  originally  at  0  now  equals  .5,  the  same 
as  the  eigenvalue  in  the  lower  right  corner.  Thus,  we  can  no  longer  say  that  this  matrix 
has  a  cluster  of  eigenvalues  near  0  and  one  near  .5,  and  so  we  can  no  longer  talk  about  the 
sensitivity  of  the  mean  of  the  cluster.  We  can  also  no  longer  identify  a  unique  eigenvector 


associated  with  an  eigenvalue  near  .5;  the  eigenvectors  have  become  ill-posed.  Indeed,  with 
additional  arbitrarily  small  perturbations  the  two  eigenvectors  for  the  eigenvalues  at  .5  can 
be  made  to  rotate  arbitrarily  within  a  two  dimensional  subspace,  or  one  of  them  can  even 
disappear.  Thus,  only  if  we  bound  \\E\\  to  be  some  value  less  than  2~^°  can  we  hope  to 
have  error  bounds.  For  this  example,  the  upper  bound  computed  by  our  software  will  be 
approximately  2  ■  10"''  (Bound  1).  For  ||£||  <  2  •  10~^,  our  upper  bound  on  the  change  in 
the  mean  of  the  eigenvalue  cluster  will  be  2l|£'||  (Bound  5),  and  our  bound  on  how  much  the 
right  invariant  subspace  can  rotate  will  be  arctan(2731||£||/(l  -  5462||i^||))  radians  (Bound 
7),  both  close  to  the  asymptotic  bounds. 

In  this  example,  it  is  easy  to  identify  the  clusters  by  inspection  of  Aq.  This  is  not 
always  the  case  in  practice,  when  the  user  is  confronted  with  a  matrix  whose  eigenvalues 
form  a  cloud  rather  than  well  separated  clusters.  Unfortunately,  there  is  as  yet  no  reliable, 
automated  procedure  for  clustering  eigenvalues;  see  [8,9,28]  for  discussion.  Our  software 
merely  provides  the  tools  for  evaluating  a  particular  clustering.  A  good  cluster  will  have  a 
much  less  sensitive  mean  and  invariant  subspace  than  any  subcluster,  and  must  be  made 
part  of  a  much  larger  (or  trivial)  cluster  before  it  becomes  significantly  less  sensitive.  The 
10  zero  eigenvalues  of  Aq  satisfy  this  criterion. 

There  is  a  very  large  literature  on  perturbation  theory  for  the  eigenproblem.  See 
[3,8,9,11,15,16,22,24,25,27,28]  for  various  theoretical  bounds.  Chan,  Feldman  and  Parlett[6] 
provided  a  Fortran  routine  to  compute  the  condition  number  of  simple  eigenvalue  in  con- 
junction with  EISPACK  routines  ORTHES  and  HQR,  but  it  does  not  provide  any  information 
about  conditioning  for  eigenvectors  and  subspaces.  Ruhe[20]  suggested  using  the  Golub- 
Reinsch  SVD  algorithm  to  calculate  the  condition  number  for  eigenvectors,  but  this  requires 
0{n^^)  flops  per  eigenvalue-eigenvector  pair,  which  is  too  expensive.  Van  Loan[26]  developed 
an  efficient  algorithm  for  estimating  condition  numbers  of  all  eigenvalue-eigenvector  pairs  of 
a  Hessenberg  matrix.  It  only  costs  0{n'^)  flops  per  eigenpair,  assuming  that  the  eigenvalues 
are  known. 

We  have  developed  new  algorithms,  which  assume  the  matrix  has  been  reduced  to  Schur 
canonical  form  (real  or  complex).  Reduction  to  Schur  canonical  form  is  done  by  LAPACK 
subroutines  SGEHRD  and  SHSEQR  in  the  real  case,  and  CGEHRD  and  CHSEQR  in  the  complex 
case.  Since  this  reduction  is  done  via  orthogonal  (or  unitary)  similarities,  the  condition 
numbers  are  identical  to  those  of  the  original  matrix.  As  we  will  see,  starting  with  the 
matrix  in  Schur  form  simplifies  many  of  the  algorithms  and  lets  us  use  existing  condition 
estimation  software  for  (quasi)triangular  matrices  [14,18,19]. 

The  rest  of  this  report  is  organized  as  follows.  Section  2  discusses  spectral  projectors  and 
the  separation  of  matrices,  quantities  on  which  later  bounds  are  based.  Section  3  discusses 
the  upper  bound  on  ||-£||f  for  our  global  error  bounds.  Section  4  discusses  asymptotic  and 
global  bounds  for  eigenvalues  and  means  of  clusters  of  eigenvalues.  In  section  5,  we  first 
define  the  angle  between  two  subspaces,  the  quantity  bounded  by  our  error  bounds.  Second, 
we  present  asymptotic  and  global  perturbation  bounds  for  both  right  eigenvectors  and  right 
invariant  subspaces.  Third,  we  discuss  (block)diagonaljzing  a  matrix  by  a  similarity.  The 
results  in  sections  2  through  5  are  stated  without  proof;  references  to  proofs  in  the  literature 
are  given.  A  tabular  summary  of  all  bounds  is  given  in  section  6.  Sections  7  and  8  describe 
the  usage  of  the  LAPACK  routines  STRSNA  and  STRSEN  for  estimating  the  desired  condition 
numbers  (actually  their  reciprocals).    STRSNA  computes  the  reciprocal  condition  numbers 


of  user-specified  eigenvalues  and/or  eigenvectors  of  tlie  input  matrix.  STRSEN  computes 
the  reciprocal  condition  numbers  of  the  mean  and/or  invariant  subspace  of  a  single  user- 
specified  cluster  of  eigenvalues.  Two  examples  are  provided  to  show  how  to  use  these  codes. 
Outlines  of  the  algorithms  are  also  given. 

The  first  two  appendices  describe  details  of  the  solution  of  the  Sylvester  matrix  equation 
and  swapping  diagonal  blocks  of  a  quasitriangular  matrix.  The  third  appendix  lists  the 
names  and  basic  functions  of  LAPACK  routines  needed  for  tlie  nonsymmetric  eigenvalue 
problem. 

We  end  with  some  notation  we  will  need  later.  Capital  letters  are  used  to  denote 
matrices,  the  corresponding  lowercase  letter  with  the  subscript  ij  referring  to  the  (i.j) 
component  (e.g.,  a,j  is  the  (z,  j)  component  of  ,4).  A  submatrix  of  a  matrix  A  is  written  as 
Aij.  Vectors  are  aJso  denoted  by  lowercase  letters  and  will  be  clearly  indicated  in  the  text. 
Lowercase  Greek  letters  will  denote  scalars. 

||x||i,  ||x||2  and  l|2:||co  denote  the  one-norm,  the  Euclidean  norm,  and  the  infinity-norm, 
respectively,  of  the  n-vector  x: 

n  /   n  \  1/2 

II     11  V^ I  /  V^  I      i2  \  II     II  II 

IfIIi  =  Z^F>K      IfI|2=     Z^F.i  ,      ||x||oo  =   max  |x,|. 

1=1  \t=i         / 

ll^lli)  ll^lb)  Halloo,  ||j^||f  denote  the  matrix  norms: 


l<:<n 


|r||i  =  max^|/,j|,      ||r||2  =  sup 


\Tx\ 


x^O 


||T||oo  =  inax^|i.,|,      IITIIf 

Note  that  ||    ||2  and  ||    1|f  are  invariant  with  respect  to  unitary  transformation. 

We  will  throughout  let  €2  denote  ||.E||2,  and  ep  denote  H^Ufi  the  norms  of  our  pertur- 
bation matrix. 

The  condition  number  of  T  is  k{T)  =  ||r||2||r~^||2.  A  subspace  spanned  by  the  columns 
of  matrix  A  is  denoted  as  R{A)  (the  range  of  matrix  A).  X{A)  denotes  the  set  of  all 
eigenvalues  of  matrix  A.  A  ®  B  denotes  the  Kronecker  product  of  two  matrices:  A  ®  B  = 
{a„B). 

The  Schur  matrix  (or  Schur  form)  of  a  real  matrix  is  an  orthogonally  similar  upper 
quasi-triangular  matrix  whose  2  by  2  diagonal  blocks  (if  any  exist)  are  of  the  form 

a    13 

Such  a  block  has  complex  conjugate  eigenvalues  q  ±/i  where  fi'^  =  -Pf.  The  Schur  form  of 
a  complex  matrix  is  a  unitarily  similar  upper  triangular  matrix. 

2      Spectral  Projectors  and  the  Separation  of  Two  Matrices 

To  explain  the  bounds  in  later  sections,  we  need  to  introduce  two  quantities,  the  spectral 
projector  P  [22,15],  and  the  separation  of  two  matrices  A  and  B,  sep(i4,i?)  [22]. 


Suppose  our  cluster  contains  m  >  1  eigenvalues,  counting  multiplicities.  Assume  the  n 
by  71  matrix  .4  is  in  Schur  canonical  form 


0       A22 


(2.1; 


where,  the  eigenvalues  of  the  m  by  m  matrix  .4ii  are  exactly  those  in  which  we  are  interested. 
In  practice,  if  the  eigenvalues  on  the  diagonal  of  A  are  in  the  wrong  order,  they  are  sorted 
to  put  the  desired  ones  in  the  upper  left  corner  as  shown  by  using  the  subroutine  STREXC 
in  Appendix  B. 

We  define  the  spectral  projector,  or  simply  projector  P  belonging  to  the  eigenvalues  of 
All  ^s 

Im       R 

0       0 


P  = 


(2.2) 


where  R  satisfies  the  system  of  linear  equations 


AuR  -  RA22  -  A12 


(2.3) 


Equation  (2.3)  is  called  a  Sylvester  equation.  Given  the  Schur  canonical  form  (2.1),  we 
solve  the  Sylvester  (2.3)  for  R  using  subroutine  STRSYL  in  Appendix  A. 

P  has  several  important  properties.  First,  the  space  spanned  by  its  columns  is  the  same 
as  the  right  invariant  subspace  of  A  belonging  to  An.  Second,  the  space  spanned  by  its 
rows  is  the  same  as  the  left  invariant  subspace  of  A  belonging  to  A\\.  Thus,  P  describes  the 
spaces  spanned  by  both  the  left  and  right  eigenvectors  belonging  to  A\\.  Third,  its  norm 
||-f||2  =  (1  +  ll-R||2)^      plays  an  important  role  in  our  error  bounds,  as  we  will  see. 

In  practice,  we  do  not  use  ||P||2  for  m  >  1  since  this  is  expensive  to  compute,  but  rather 
the  cheaper  overestimate 


\P\\'^{l^\\R\\lYl^ 


(2.4) 


The  separation  sep(Aii,yl22)  of  the  matrices  An   and  A22  is  defined  as  the  smallest 
singular  value  of  the  linear  map  in  (2.3)  which  takes  A'  to  A\\X  —  X A22,  i-e. 


Sep(Aii,A22) 


mm 


l^iiA  -  A'A22||f 
liA'llF 


(2.5) 


This  formulation  lets  us  estimate  sep  using  the  condition  estimator  SLACON  [14,18,19],  which 
estimates  the  norm  of  a  linear  operator  HTHi  given  the  ability  to  compute  Tx  and  T^y 
quickly  for  arbitrary  x  and  y.  In  our  case,  multiplying  an  arbitrary  vector  by  T  means 
solving  the  Sylvester  equation  (2.3)  with  an  arbitrary  right  hand  side,  and  multiplying  by 
T-^  means  solving  the  same  equation  with  A\i  replaced  by  AJ^  and  A22  replaced  by  A2'2- 
Solving  either  equation  costs  at  most  0{rfi)  operations,  or  as  few  as  0{n^)  if  m  <C  n. 
Another  formulation  which  in  principle  permits  an  exact  evaluation  of  sep(Aii,  A22)  is 


Sep(Aii,  A22)  =  <7inin(/n-m  ®  ^11  "  ^22  ®  Im) 


(2.6) 


where  o^in  is  the  smallest  singular  value.   This  method  is  generally  impractical,  however, 
because  the  matrix  whose  smallest  singular  value  we  need  is  m{n  —  m)  dimensional,  which 


can  be  as  large  as  7?-/4.  Thus  we  would  require  as  much  as  0{n'^)  extra  workspace  and 
0{n^)  operations,  much  more  than  the  estimation  method  of  the  last  paragraph. 

sep(Aii,  A22)  measures  the  "separation"  of  the  spectrum  of  An  and  A22  in  the  following 
sense.  It  is  zero  if  and  only  if  An  and  A22  have  a  common  eigenvalue,  and  small  if  there  is 
a  small  perturbation  of  either  one  that  makes  them  have  a  common  eigenvalue.  If  An  and 
,422  are  both  normal  matrices,  then  se-p(Au,A22)  is  just  the  minimum  distance  between  an 
eigenvalue  of  Au  and  an  eigenvalue  of  A22- 

STRSNA  computes  1/||P||2  (which  is  always  <  1,  avoiding  the  possibility  of  overflow)  and 
sep  for  user-selected  individual  eigenvalues  (i.e.  Au  is  1  by  1).  STRSEN  computes  l/||P|r 
and  sep  for  a  single  user-specified  cluster  with  m  >  1  eigenvalues. 

3  An  Upper  Bound  on  \\E\\f  for  Global  Error  Bounds 

We  discuss  the  upper  bound  on  \\E\\f  which  limits  the  applicability  of  our  global  bounds 
in  the  next  two  sections.  As  stated  in  the  introduction,  this  upper  bound  occurs  because  if 
||i^||f  is  large  enough  that  the  eigenvalue  being  considered  (or  one  of  the  eigenvalues  in  the 
cluster  being  considered)  moves  and  coalesces  with  another  eigenvalue  (outside  the  cluster), 
then  we  can  no  longer  uniquely  identify  the  cluster  for  which  we  want  bounds.  Thus,  in 
this  section  we  present  lower  bounds  on  the  smallest  ||£^||f  such  that  A  +  E  has  a  multiple 
eigenvcdue  (or  a  multiple  eigenvalue  involving  at  least  one  eigenvalue  within  the  original 
cluster  and  one  outside). 

Bound  1:  [22,  Theorem  4.14]  Let  A,  P  and  sep{An,A22)  be  defined  as  in  (2.1),  (2.2)  and 
(2.5).   Then  as  long  as 

,  sep(A„,A22)  ,      . 

I'^"^<       4.IIPII2  ^^-'^ 

the  eigenvalues  in  the  cluster  belonging  to  An  will  remain  disjoint  from  the  eigenvalues  out- 
side the  cluster.  In  particular,  the  global  error  bounds  of  sections  4  o-nd  5  will  be  guaranteed 
valid  only  for  E  satisfying  (3.7).  We  may  replace  \\P\\2  by  \\P\\'  as  defined  in  (2-4)  to  get 
a  slightly  smaller  upper  bound. 

Bound  1  can  be  quite  conservative,  greatly  underestimating  the  smallest  perturbation 
needed  to  make  eigenvalues  coalesce.  However,  it  is  nearly  exact  in  some  cases  (e.g.  for 
2  by  2  matrices  and  normal  matrices),  and  a  good  estimate  in  many  others;  see  [8,9]  for 
discussion. 

1/||P||'  (or  I/IIPII2  if  m  =  1)  and  sep(ylii,  A22)  are  computed  by  STRSNA  and  STRSEN  as 
described  in  section  2. 

4  Conditioning  of  Eigenvalues 

In  this  section,  we  review  how  to  measure  the  sensitivity  of  both  simple  eigenvalues  and 
clusters  of  eigenvalues. 


4.1      Conditioning  of  Simple  Eigenvalues 

Let  A  be  a  simple  eigenvaJue  of  the  n  by  n  matrix  A,  with  unit  right  eigenvector  x  and  unit 
left  eigenvector  y.  In  other  words  Ax  —  Xx,  y^ A  =  Ay-^,  and  ||a:||2  =  ||2/||2  =  1-  Let  P  be 
the  spectral  projector  for  A;  one  may  write  P  =  {x  ■  y'^)/{y'^  ■  x).  Note  that  ||P||2  =  l/l2/"^^K 
the  secant  of  the  angle  between  x  and  y. 

Let  £■  be  a  perturbation  of  A,  and  fj  =  Il-^ll2-   Let  A'  be  the  perturbed  eigenvalue  of 
A  +  E. 

Bound  2:  [27,  p.  68] 

|A'-A|   <  6211^112  +  0(6^) 


The  Oiel)  term  indicates  this  is  an  asymptotic  bound,  applicable  only  for  sufficiently 
small  €2-  This  bound  is  attainable,  in  the  sense  that  for  €2  sufficiently  small,  there  exists 
an  E  such  that  |A'  -  A|  =  e2||P||2  +  0(62)- 

There  is  also  a  global  version  of  this  bound: 

Bound  3:  [3]  Suppose  A  has  all  simple  eigenvalues  A,.  Let  P,  be  the  corresponding  spectral 
projectors.   Then  any  eigenvalue  A'  of  A  +  E  must  lie  in  one  of  the  disks 

{A:     |A-A.|<n€2||P.||2} 


There  is  no  limit  on  the  size  of  €2  in  Bound  3.  Note  that  the  sizes  of  the  disks  are  just 
n  times  larger  than  in  Bound  1,  where  €2  must  be  small.  Bound  3  is  an  stronger  version  of 
what  is  often  called  the  Bauer-Fike  Theorem,  although  Bauer  and  Fike  proved  this  stronger 
version  as  well.  In  the  weaker  Bauer-Fike  Theorem,  all  of  the  disks  have  the  same  radius, 
approximately  equal  to  the  largest  radius  max,  nc2||P:||2  in  Bound  3.  Note  that  Bound  3 
is  only  useful  when  all  the  radii  ne2||Pi||2  a^re  of  modest  size,  since  if  one  or  more  disks  is 
so  large  that  it  intersects  all  the  other  disks,  there  is  little  information  about  locations  of 
individual  eigenvalues;  we  only  know  they  lie  in  the  union  of  all  the  disks. 

The  subroutine  STRSNA  can  compute  1/||P||2  for  a  user-specified  subset  of  the  eigenvalues 
of  A. 

4.2      Conditioning  of  Clustered  Eigenvalues 

It  is  easiest  to  think  of  A  in  Schur  form  (2.1),  with  the  eigenvalues  of  An  being  the  cluster 
of  interest.  We  are  interested  in  bounding  the  perturbation  in  the  average  of  the  eigenvalues 
of  the  cluster,  which  may  be  written  tvAu/m,  the  trace  of  An  divided  by  m.  Let  E  be 
a  perturbation  of  A,  and  €2  =  ll-E^lb-  Let  A  =  tr/ln/m  be  the  mean  of  the  unperturbed 
eigenvalues,  and  A'  be  the  mean  of  the  perturbed  eigenvalues. 

Bound  4:  [15,  sec.  II. 2.2] 

jA-  A'j   <  C2||P||2  +  0(C2) 


We  may  substitute  \\P\\'  of  equation  (2-4)  for  \\P\\2  to  get  a  slightly  weaker  bound. 

The  Oiel)  indicates  tliis  is  an  asymptotic  bound,  applicable  only  for  sufficiently  small 
€2-  This  bound  is  nearly  attainable,  in  the  sense  that  for  (2  sufficiently  small,  there  exists 
an  E  such  that  |A  -  A'|  <  ||P||2£2/"^  +  ^(ej).  When  m  =  1,  it  is  of  course  identical  to  the 
bound  in  the  previous  subsection. 

Our  global  bound  on  |A  —  A'|  will  only  be  valid  for  H-EHf  satisfying  Bound  1  of  section 
2: 

Bound  5:  [22,  page  748]  Suppose  \\E\\p  satisfies  Bound  1.   Then 

|A-  A'l  <  2e2|lP||2 

Thus,  the  global  bound  is  just  twice  as  large  as  the  asymptotic  bound.  Again  we  may  sub- 
stitute \\F\\'  of  equation  (2.4)  for  \\P\\2  to  get  a  slightly  weaker  bound. 

STRSNA  computes  1/||P||2  for  a  user-specified  set  of  individual  eigenvalues.  STRSEN  can 
compute  1/||P||'  for  a  single  user-specified  cluster  of  m  >  1  eigenvalues. 

4.3      Stability 

When  the  eigenvalues  of  a  full  matrix  A  have  been  found  from  its  computed  Schur  form  A, 
the  computed  s  will  be  those  appropriate  to  A.  These  s  can  differ  substantially  from  the 
true  s.  Indeed,  when  A  is  defective,  A  will  usually  not  be,  and  hence  zero  s  will  become 
nonzero  s.  The  reverse  situation  could  occur  though  this  is  much  less  probable.  Since  some 
of  the  s  may  not  even  be  the  correct  order  of  magnitude,  it  might  be  felt  that  our  heavy 
reliance  on  them  is  unjustified.  Since  our  algorithms  for  computing  the  Schur  form  and  5(A) 
are  backward  stable,  s{\)  is  the  correct  value  for  a  matrix  A  ■\-  E  very  close  to  the  original 
matrix  A.  This  justifies  their  use.  The  same  comments  apply  to  the  computation  and  use 
of  sep  described  in  the  next  section. 

5      Conditioning  of  Right  Eigenvectors  and  Right  Invariant 
Subspaces 

In  this  section,  we  review  how  to  measure  the  sensitivity  of  eigenvectors  and  invariant 
subspaces.  We  begin  by  defining  of  the  angle  between  two  subspaces,  and  then  use  it  to 
describe  the  conditioning  of  eigenvectors  and  invariant  subspaces. 

5.1      Angles  Between  Subspaces 

Let  0{x,y)  denote  the  acute  angle  between  two  nonzero  n- vectors  x  and  y.  We  may  write 
cos6{x,y)  =  \x'^y\/{\\x\\2\\y\\2)-  We  wish  to  generalize  this  to  the  (maximum)  angle  between 
two  TO  >  1  dimensional  subspaces,  which  we  denote  A'  and  y.  One  way  to  define  this  angle 
is  as 

^max(A',;y)  =     max        min     9(x,y)  {-     max        min     d{x,y))  (5.8) 

X  e  A'    y  £  y  y  ey    X  e  A' 

x^O      2//0  J//0      I/O 
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A  more  computational  definition  of  ^maxl'^'i^)  is  the  following.  Suppose  A'  is  spanned 
by  the  columns  of  the  n  by  m  orthonormal  matrix  A',  and  similarly  y  is  spanned  by  the 
columns  of  the  n  by  m  orthonormal  matrix  Y .  Then 

einax(A',>')  =  arccosan^n(r^A')  (5.9) 

Our  bounds  of  the  next  two  sections  will  bound  Onieix{A\A")  where  A'  is  an  unperturbed 
invariant  subspace,  and  A"  is  a  perturbed  invariant  subspace. 
We  may  also  define  the  minimum  angle  between  A'  and  y  as 

^min('^'i3-')  =      min         min     0{x,y)   (=      min         min      6(x,y)) 

X  e  A'    y  ey  y  ey    X  e  A' 

X  ^0      y  j!:0  y  :fLO      X  yiO 

This  may  be  reexpressed  computationally  as 

^Tnini^^y)  -  arCCOSCrmax(5''^A') 

The  norms  of  the  spectral  projectors  ||P||2  introduced  in  section  2  have  a  simple  inter- 
pretation in  terms  of  angles  between  subspaces.  Let  P  be  the  spectral  projector  for  the 
eigenvalue  cluster  with  right  invariant  subspace  7v  and  left  invariant  subspace  C  Let  TZc  be 
the  complementary  right  invariant  subspace  (the  subspace  for  the  other  eigenvalues)  and 
Cc  be  the  complementary  left  invariant  subspace.  Then 

II-PII2   =   csce^jn,n,)  =  CSC  e^nicXc) 

\\P\\2     =     sec^™^^(7^,£)  =  sece^.ax(7^c,-Cc) 

In  other  words,  as  ||P||2  grows  and  the  cluster  becomes  more  ill-conditioned,  the  minimum 
angle  between  complementary  right  (or  complementary  left)  subspaces  shrinks.  Also,  the 
maximum  angle  between  corresponding  left  and  rigth  invariant  subspaces  grows. 

5.2      Conditioning  of  Right  Eigenvectors  and  Right  Invariant  Subspaces 

We  assume  A  is  in  Schur  canonical  form  (2.1),  with  the  eigenvalues  of  An  being  the  cluster 
whose  right  invariant  subspace  TZ  we  are  interested  in.  Let  £  be  a  perturbation  of  A,  and 
ff  =  II^IIf-  Let  7v^'  be  the  perturbed  right  invariant  subspace  of  A  +  E. 


Bound  6:  [8,  Lemma  7.8] 


^max(7^,  7^')  <  — ^  f  ^  ^      ^  +  0(4) 


The  O(e^)  indicates  that  this  is  an  asymptotic  bound,  applicable  only  for  sufHciently 
small  c/r.  It  is  nearly  attainable  for  sufficiently  small  ep. 

Bound  7:  [8,  Lemma  7.8]  Suppose  \\E\\f  satisfies  Bound  1.   Then 

^max(^5^')  <  arctan  f 


sep(/lii,A22)-4||P||2eF 
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Il-Plb  ^cy  be  replaced  by  \\P\\'  of  equation  (2.4)  to  obtain  a  slightly  weaker  bound. 

Bounds  6  and  7  imply  that  sep  is  the  reciprocal  of  the  condition  number  for  eigen- 
vectors and  invariant  subspaces.  Routines  STRSNA  and  STRSEN  compute  sep  for  individual 
eigenvectors  and  a  given  invariant  subspace,  respectively. 

5.3      (Block)diagonalizing  a  Matrix  with  a  Similarity 

Occasionally  one  wishes  to  find  a  matrix  V  which  diagonalizes  A:  V~^AV  =  A,  where  A  is 
a  diagonal  matrix  with  the  eigenvalues  on  the  diagonal.  This  may  be  useful  for  computing 
functions  of  matrices.  For  example,  to  exponentiate  a  matrix  one  may  use  the  identity 
exp(i4)  =  V'exp(A)y~''  =  Vdiag{e^^, . .  .,e'^")V~^ .  The  accuracy  of  such  an  algorithm 
depends  on  the  condition  number  k{V)  of  V .  The  columns  of  V  must  be  eigenvectors  of  A, 
but  their  norms  are  arbitrary;  we  would  like  to  choose  these  norms  to  minimize  k(T').  The 
next  bound  gives  a  nearly  optimal  choice  of  the  norms  of  the  columns  of  V ,  and  bounds  the 
resulting  ii{V). 

Bound  8:  [7]  Suppose  A  has  distinct  eigenvalues  A,  with  corresponding  right  eigenvectors 
V,,  where  we  assume  ||v,||2  =  1,  and  projectors  P,.  Let  V  =  [ui,...,v„]  be  the  matrix  of 
these  eigenvectors.  Let  V  =  [qjUi,  . . . ,  a„t'n]  be  another  matrix  where  the  a;  are  arbitrary 
nonzero  scalars.   Then 

inax||P,||2  <  k(V') 

i 

Also 

max||P,||2  <  k{V)  <  n-max||P,||2 

1  : 

In  other  words  choosing  the  columns  of  V  to  have  norm  1  nearly  minimizes  k{V)  over  all 
matrices  whose  columns  are  eigenvectors. 

A  variation  on  diagonalizationis  block-diagonalization,  where  we  ask  only  that  V~^  AV  = 
B  be  block  diagonal,  where  the  diagonal  blocks  P,,  of  D  have  specified  eigenvalues  (which 
are  all  disjoint  subsets  of  the  eigenvalues  of  A).  Suppose  P„  is  located  in  rows  and  columns 
j  through  k  of  B.  Then  columns  j  through  k  of  V  must  span  the  right  invariant  subspace 
of  A  corresponding  to  the  eigenvalues  of  B.  Let  V;  denote  these  columns  of  V .  Just  as  we 
could  choose  the  norm  of  each  column  of  V  when  we  wanted  to  diagonalize  A,  here  we  have 
the  freedom  to  choose  F,  to  be  any  basis  of  the  right  invariant  subspace  we  like.  Again,  we 
would  like  to  choose  the  basis  which  minimizes  k{V).  The  next  bound  says  how  to  do  this. 

Bound  9:  [7]  Let  the  set  \{A)  of  eigenvalues  of  A  be  written  as  a  disjoint  union  uj'_j5,  ofb 
sets  of  eigenvalues  Si.  Let  n,  be  the  number  of  eigenvalues  in  S,,  counting  multiplicities.  Let 
Pi  be  the  projector  corresponding  to  S,,  and  7v,  the  corresponding  right  invariant  subspace. 
Let  Vi  be  any  matrix  whose  Ui  columns  form  an  orthonormal  basis  of  7Z,,  and  write  V  — 
[Vi, . . .,  V(,].  Then  V~^AV  =  B  is  block  diagonal  where  diagonal  block  P,,  has  eigenvalues 
Si.  Let  V-'  be  an  arbitrary  matrix  whose  n,  columns  form  a  basis  o/7v,,  and  write  V  = 
[Vj', . . . ,  Vj,'].  V'~^AV'  =  B'  is  also  block  diagonal  where  diagonal  block  P,',  has  eigenvalues 
5,.   Then 

max||P,||2  <  H.{V') 
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Also 


max  \\Pi\\2  <  k(  V)  <  b  ■  max  ||P,||2 


In  other  words  choosing  the  columns  o/K  which  span  T^i  to  be  orthonormal  nearly  minimizes 
k(V)  over  all  block-diagonalizing  similarities  which  reduce  A  to  diagonal  blocks  with  the  same 
eigenvalues  in  each  block. 

6      Summary:   Perturbation  Table 

For  convenience,  the  bounds  presented  in  tlie  preceding  sections  are  summarized  in  the 
following  table.  The  notation  is  as  follows.  We  assume  the  matrix  A  is  in  Schur  canonical 
form  (2.1).  P  denotes  the  spectral  projector  associated  with  with  eigenvalue(s)  of  An 
defined  in  (2.2).  sep  will  be  shorthand  for  sep(^ii,  >l22),  defined  in  (2.5).  A  will  denote  the 
unperturbed  eigenvalue  if  A\\  is  1  by  1,  and  if  An  is  larger  A  will  denote  the  unperturbed 
mean  of  its  eigenvalues.  A'  and  A'  will  denote  the  perturbed  vaJues  of  A  and  A,  respectively, 
for  A-'r  E.  TZ  denotes  the  unperturbed  right  invariant  subspace  corresponding  to  An,  and 
7v'  denotes  its  perturbed  counterpart  oi  A  ■\-  E.  6  will  denote  ^max(^i7v'),  the  angular 
perturbation  of  the  right  invariant  subspace  as  defined  in  (5.8)  or  (5.9).  £2  will  denote 
||£'||2  and  ep  will  denote  II^^Ufi  norms  of  the  perturbation  E.  In  the  table,  each  asymptotic 
bound  has  a  +0(e^)  term  which  is  not  written.  Superscripts  in  parentheses  on  each  bound 
indicate  which  Bound  in  the  body  of  text  they  are.  The  superscript  f  indicates  that  the 
bound  applies  only  li  ep  <  sep/(4||P||2)  (Bound  1). 


Asymptotic  Bounds 

Global  Bounds 

Simple  eigenvalue 

|A-A'|<e2||P||2(2) 

|A-A'|<ne2||P||2'" 

Eigenvalue  Cluster 

|A-A'|<.2||P||2*'*^ 

|A-A'|<2e2i|P||2^<'' 

Invariant  subspace 

"  -    sep 

^<arctan(,^4,P|,^,J^" 

In  addition,  Bound  8  says  that  the  nearly  best  conditioned  matrix  V  such  that  V~^AV 
is  diagonal  has  as  its  columns  the  eigenvectors  of  A  all  with  norm  equal  to  1.  The  condition 
number  n{V)  of  this  V  satisfies  max,  ||P,||2  <  k,{V)  >  n  ■  max,  ||P,||2,  where  P,  is  the 
projector  corresponding  to  eigenvalue  A,. 

Bound  9  describes  a  nearly  best  conditioned  matrix  V  such  that  V~^AV  =  B  is  block 
diagonal,  such  that  the  b  diagonal  blocks  of  B  have  specified  eigenvalues.  This  nearly 
best  V  may  be  written  V  =  [Vi,...Vb]  where  1-^  is  any  orthonormal  basis  of  the  right 
invariant  subspace  of  A  corresponding  to  the  eigenvalues  in  the  i-th  diagonal  block  of  B. 


13 


The  condition  number  k(\')  of  1'  satisfies  max,  ||Pi||2  <  k(^')  <  i^  ■  max,  ||P,||2-  wliere  P,  is 
the  projector  corresponding  to  eigenvalues  of  the  7-th  diagonal  block  of  B. 

In  summary,  we  see  that  all  the  condition  numbers  for  the  nonsymmetric  eigenproblem 
depend  on  the  quantities  ||-P||2  and  sep.  The  use  LAPACK  subroutines  STRSNA  and  STRSEN 
to  estimate  these  quantities  is  discussed  in  the  following  sections. 

7     STRSNA  -  Estimating  the  Condition  of  Individual  Eigen- 
pairs 

In  this  section,  we  show  how  to  use  LAPACK  subroutine  STRSNA  for  estimating  condition 
numbers  of  selected  eigenvalues  and  eigenvectors.  We  also  describe  the  algorithms  used.  The 
variable  SEP  corresponds  to  the  sep  of  the  preceding  sections,  and  the  variable  S  corresponds 
to  1/||P||2.  We  return  1/||P||2  instead  of  ||P||2  since  1/||P||2  is  between  zero  and  one.  thus 
avoiding  overflow  problems  for  very  ill-conditioned  eigenvalues. 

7.1      Usage 
Single  precision. 

CALL   STRSNA (SELECT , N , T , LDT , RE , LORE , LE , LDLE , S , SEP , 

$  MM. M, WORK. LDWORK.X.V.B.ISGN. INFO)    • 

* 

*  . .  Scalar  Arguments  . . 

INTEGER  N. LDT. LORE, LDLE, MM. M.LDWORK, INFO 

♦ 

*  . .  Array  Arguments  . . 
LOGICAL  SELECTC*) 
INTEGER  ISGN(*) 
REAL            S(*),  SEP(*) 

REAL  T(LDT.*).  RE(LDRE.*),  LE(LDLE.*) 

REAL  WORKCLDWORK,*)  .  X(*)  ,  V(*)  ,  B(*) 

* 

*  Arguments 

*  ========= 

* 

*  SELECT  -  LOGICAL     array  if  DIMENSION  (N) . 

*  On  entry.  SELECT  specifies  the  eigenpair  whose 

*  condition  numbers  are  to  be  estimated.  The  condition 

*  number  corresponding  to  the  J-th  eigenpair  is  specified 

*  by  setting  SELECT(J)  to  .TRUE.. 

*  Not  modified. 
* 

*  N      -  INTEGER 

*  On  entry,  N  specifies  the  order  of  the  matrix  T.  N  must 

*  be  at  least  zero. 
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* 

* 

* 

* 

* 

* 
* 

* 

* 
* 
* 
* 
* 
* 
* 
* 
* 

* 

* 
* 
* 
* 
* 

* 

* 
* 
* 

* 
* 
* 

* 


Not  modified. 

T      -  REAL  array  of  DIMENSION  (LDT.N) . 

On  entry,  T  contains  an  upper  quasi-triangular  matrix  in 
Schur  canonical  form.  This  means  that  the  diagonal  entries 
of  2  by  2  diagonal  blocks  must  be  equal. 
Not  modified. 

LDT    -  INTEGER 

On  entry,  LDT  specifies  the  first  dimension  of  T  as 
declared  in  the  calling  (sub)program.  LDT  must  be  at 
least  maxd  ,  N)  . 
Not  modified. 

RE     -  REAL  array  of  DIMENSION  (LDRE.MM). 

On  entry,  RE  contains  the  real  and  imaginary  parts  of  the 
selected  right  eigenvectors  computed  by  STREVC  or  SHSEIN. 
If  the  next  selected  eigenvalue  is  real,  the  next  col\imn 
of  RE  contains  its  eigenvector.  If  the  next  selected 
eigenvalue  is  complex,  the  next  two  coltimns  of  RE  contain 
the  real  and  imaginary  parts  of  its  eigenvector. 
Not  modified. 

LORE   -  INTEGER 

On  entry,  LDRE  specifies  the  leading  dimension  of  RE  as 
declared  in  the  calling  (sub)program.  LDRE  must  be  at 
least  maxd  ,  N)  . 
Not  modified. 

LE     -  REAL  array  of  DIMENSION  (LDLE.MM). 

On  entry,  LE  contains  the  real  and  imaginary  parts  of  the 
selected  left  eigenvectors  computed  by  STREVC  or  SHSEIN. 
If  the  next  selected  eigenvalue  is  real,  the  next  column 
of  RE  contains  its  eigenvector.  If  the  next  selected 
eigenvalue  is  complex,  the  next  two  columns  of  LE  contain 
the  real  and  imaginary  parts  of  its  eigenvector. 
Not  modified. 

LDLE   -  INTEGER 

On  entry,  LDLE  specifies  the  leading  dimension  of  LE  as 
declared  in  the  calling  (sub)program.  LDLE  must  be  at 
least  maxd  ,  N)  . 
Not  modified. 


REAL 


array  of  DIMENSION (MM) 
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*  Dn  exit,  S  contains  the  reciprocals  of  the  condition 

*  numbers  of  the  selected  eigenvalues.  If  the  Jth  and 

*  (J+l)st  eigenpairs  are  complex  conjugate,  then  both 

*  S(J)  and  S(J+1)  will  be  set  (and  equal). 
* 

*  SEP    -  REAL  array  of  DIMENSION(MM) . 

*  On  exit,  SEP  contains  the  estimated  reciprocals  of  the 

*  condition  numbers  of  the  selected  right  eigenvectors. 

*  If  the  Jth  and  (J+l)st  eigenpairs  are  complex  conjugate, 

*  then  both  SEP(J)  and  SEP(J+1)  will  be  set  (and  equal) . 
* 

*  MM     -  INTEGER 

*  Dn  entry,  MM  should  be  set  to  an  upper  bound  for  the 

*  length  of  arrays  S(*)  and  SEP(*)  required  to  store  the 

*  reciprocal  condition  numbers  to  be  estimated.  Note  that 

*  for  a  complex  conjugate  eigenpair,  we  need  two  locations 

*  for  S  and  SEP.  This  means  S(J) ,  SEP(J) ,  RE(J) ,  and  LE(J) 

*  correspond  to  the  same  eigenvalue  for  all  J. 

*  Not  modified. 
* 

*  M      -  INTEGER  , 

*  On  exit,  M  is  the  size  of  arrays  S(*)  and  SEP(*)  actually 

*  used  to  store  condition  numbers. 

* 

*  WORK   -  REAL  array  of  DIMENSION(LDWORK.N) 

*  Workspace. 

* 

*  LDWORK  -  INTEGER 

*  On  entry,  LDWORK  specifies  the  leading  dimension  of  WORK 

*  as  declared  in  the  calling  (sub)program.  LDWORK  must  be 

*  at  least  max(l,  N) . 

*  Not  modified. 
* 

*  X      -REAL  array  of  DIMENSI0N(2*(N-1)) . 

*  Workspace. 
* 

*  V      -  REAL  array  of  DIMENSI0N(2*(N-1)) . 

*  Workspace . 

* 

*  B      -  REAL  array  of  DIMENSION (N) 

*  Workspace. 
* 

*  ISGN   -  INTEGER  array  of  DIMENSI0N(2*(N-1)) 

*  Workspace. 
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* 

INFO 

-  INTEGER 

* 

On  exit 

* 

0 

* 

-K 

* 

N+1 

* 

INFO   is   set  to 

for  normal  return, 

input   argument  number  K   is    illegal. 

the  assigned  length  of   S  and  SEP  too   small. 


Double  precision.  The  calling  sequence  of  the  double  precision  routine  DTRSNA  is  the 
same  as  that  of  the  corresponding  single  precision  "S"  subroutine  except  that  all  the  real 
variables  are  double  precision. 

Complex.  The  calling  sequence  of  the  single  precision  complex  is  essentially  the  same 
as  STRSNA,  except  that  the  T,  RE,  WORK,  X,  V  arrays  are  complex,  and  the  integer  array 
ISGN  is  real. 

Double  precision  complex.  The  calling  sequence  of  the  double  precision  complex 
routine  ZTRSNA  is  the  same  as  that  of  the  corresponding  single  precision  "C"  subroutine 
except  that  all  the  real  variables  are  double  precision  and  all  the  complex  variables  are 
double  precision  complex. 

7.2     Example 

The  following  program  segment  illustrates  the  use  of  the  single  precision  subroutine  to 
estimate  selected  reciprocal  condition  numbers  of  the  eigenvalues  and  eigenvectors  of  a 
general  matrix.  The  program  first  reduces  a  general  matrix  to  upper  Hessenberg  form 
by  SGEHRD,  and  then  computes  the  Schur  decomposition  by  the  multishift  QR  iteration 
(SHSEQR).  After  that  the  user  should  input  the  logical  array  SELECT  to  specify  the  eigenpairs 
whose  condition  numbers  will  be  estimated,  and  STREVC  is  called  to  compute  the  selected 
right  and  left  eigenvectors  and  compactly  store  them  in  array  RE  and  LE.  Finally  STRSNA  is 
called  to  return  the  desired  reciprocal  condition  numbers. 

PROGRAM   TEST 
INTEGER   ISGN (100) 
LOGICAL   SELECT(50) 

REAL      A(50,50),  WR(50) ,  WI(50),  RE(50.50),  LE(S0,50) 
REAL      S(50),  SEP(50),  RWORK(50) 

REAL   B(50),  X(IOO),  V(IOO),  WORK(50,150),  U(50,150) 
INTEGER   N,  LDA,  LORE,  LDLE.  LDWORK,  M,  I,  J,  MAXN ,  INFO 
INTEGER   NBLCKl,  NSHIFT,  NBLCK2 
REAL      DUMMY 

PARAMETER  (LDA  =  50,  LORE  =  50,  LDLE  =  50,  LDU  =  50) 
PARAMETER  (LDWORK  =  50,  MAXN  =  150) 
* 

*  Input  data: 

*  N:  the  order  of  matrix  A. 

*  A:  N  by  N  array  to  store  the  input  matrix. 

*  NBLCKl:  the  blocksize  used  in  Hessenberg  reduction. 

*  NSHIFT:  the  number  of  shifts  used  in  multishift  QR  algorithm. 
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*  NBLCK2:  the  blocksize  used  in  bulge  chasing  in  QR  iteration. 
* 

READ(*,*)  N 

DO  10  I  =  l.N 

READ(*,*)  (A(I,J),J  =  l.N) 
10      CONTINUE 

READC*.*)  NBLCKl 

READ(*,*)  NSHIFT 

READC*,*)  NBLCK2 
* 

*  Compute  Schur  decomposition. 
* 

CALL  XENVIRC 'BLOCK' .NBLCKl) 

CALL  SGEHRDCH',  N,  A,  LDA ,  DUMMY.  DUMMY.  WR,  WORK.  LDWORK. 
$  MAXN.  INFO) 

CALL  XENVIRC 'SHIFT' .  NSHIFT) 
CALL  XENVIRC 'BLOCK' .  NBLCK2) 

CALL  SHSEQRC'S'.  N.  A,  LDA.  DUMMY,  DUMMY.  WR,  WI .  U,  LDU. 
$  MAXII,  WORK,  LDWORK.  MAXN.  INFO) 

DO  20  I  =  1,  N 

WRITEC*,*)  WRCD.WICI)  .; 

20      CONTINUE 
* 

*  Input  logical  array  SELECT  to  specify  the  eigenpairs  whose 

*  condition  numbers  will  be  estimated. 

* 

DO  30  I  =  l.N 

READC*,*)  SELECTCI) 
30      CONTINUE 

* 

*  Compute  the  selected  eigenvectors 

* 

CALL  STREVCC'B',  SELECT,  N,  A,  LDA.  RE,  LDRE,  LE,  LDLE, 
$  N.  M.  RWORK.  INFO) 

* 

*  Estimate  the  selected  condition  numbers  of  eigenpairs 
♦ 

CALL  STRSNACSELECT,  N,  A,  LDA.  RE,  LDRE,  LE,  LDLE,  S,  SEP. 
$  N,  M,  WORK,  LDWORK,  X,  V,  B,  ISGN,  INFO) 

* 

*  Print  output 
* 

DO  40  I  =  l.M 

WRITEC*,*)    SCD.SEPCI) 
40  CONTINUE 
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STOP 
END 

7.3      Outline  of  the  Algorithm 

The  STRSNA  routine  is  designed  to  estimate  tlie  reciprocals  of  condition  numbers  of  the 
selected  eigenvalue-eigenvector  pairs  of  a  Schur  canonical  matrix  T. 

Logical  arrary  SELECT  specifies  the  condition  numbers  to  be  estimated.  The  arrays 
RE  and  LE  are  used  in  STREVC  to  compactly  store  the  selected  right  and  left  eigenvectors 
respectively.  Then  RE  and  LE  are  used  in  STRSNA  to  compute  the  reciprocal  condition 
number  for  individual  eigenvalues: 


where  r  and  /  are  the  right  and  left  eigenvectors  of  T  corresponding  to  the  eigenvalue  A. 
Note  that  for  complex  eigenvalues,  the  next  two  columns  of  RE  (LE)  store  the  real  and 
imaginary  parts  of  the  right  (left)  eigenvectors,  respectively.  We  see  that  computing  the 
reciprocals  of  condition  number  of  an  eigenvalue  costs  0{n)  operations. 

Variable  S(*)  contains  the  reciprocals  of  condition  numbers  of  the  selected  eigenvalues. 

For  estimating  the  reciprocals  of  condition  numbers  of  the  associated  right  eigenvectors, 
STRSNA  first  calls  subroutine  STREXC  to  swap  the  diagonal  blocks  of  matrix  T  by  orthogonal 
transformation  to  the  form: 

«^«'  =  ( ^!>'  t ) 

where  the  rii  by  ni  matrix  Tn  is  1  by  1  or  2  by  2  depending  on  whether  the  eigenvalue  is 
real  or  complex.  If  Tn  is  a  1  by  1  block  A,  we  have 

sep(A)=    min   ||(A/ -  T22)a;||2 
Ikll2=i 

If  Til  is  a  2  by  2  block,  then  we  use  a  unitary  rotation  to  triangularize  the  2  by  2  block  to 
get 

'   A      <i2 
0      f22 


whence 


sep(A)  =    min    ||(A  -  T22)A\2- 
Ikll2=i 


In  both  cases  sep  can  be  estimated  by  estimating  the  one-norm  of 

A'-i  =  (T'  -  A)-^ 
because  of  the  relationship 

^  ■|A-l||l   <   ^TVT  =   l|A'-'ll2  <   V^r^T||A'-^||i 


\/n  -  1  sep(A) 
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Estimating  the  one-norm  of  ||A'""^||i  can  been  done  by  calling  SLACON  via  a  reverse  com- 
munication interface.  This  means  one  must  provide  the  solution  vectors  x  and  y  of  the 
quasi-triangular  systems: 

A'a:  =  z,      K'^y  =  z 

where  z  is  determined  by  SLACON.  This  is  the  function  of  the  subroutine  SLAQTR.  Note  that 
A'  is  a  complex  matrix  if  A  is  a  complex  eigenvalue,  but  it  is  of  the  form 

K  =  C  +  ID 

where  the  real  part  C  is  a  real  upper  quasi- triangular  matrix,  and  the  imaginary  part  is 


D  = 


(   X      X 

X 

V 


\ 


/ 


Hence  we  can  easily  solve  the  complex  systems 

(C  +  iD)[p  +  iq)  =  (e  +  if),     (C  +  iDf{g  +  ih)  =  (e  +  if) 

in  real  arithmetic,  and  use  2(n  -  1)  length  vectors  x  =  \  and  y  =        ,     )  as  the  input 

vectors  of  SLACON.  We  aJso  use  the  fact  that  for  any  complex  matrix  C  +  iD, 


C     -D 
D      C 


li  <llC  +  ii)lli  < 


c 

D 


-D 
C 


The  cost  of  the  algorithm  depends  on  the  location  of  the  selected  eigenvalues  along  the 
diagonal  of  the  input  matrix.  Swapping  adjacent  diagonal  blocks  costs  0(n),  so  moving 
an  eigenvalue  at  diagonal  position  k  to  the  upper  left  costs  0{kn)  operations.  Since  it 
requires  about  O(n^)  to  solve  a  quasi- triangular  system,  estimating  the  condition  number 
of  an  eigenvector  costs  0(n'^)  operations  once  the  eigenvalue  is  in  the  upper  left  corner. 
Therefore  the  total  cost  is  0{n^)  per  eigenvector  condition  number. 

The  variable  SEP(*)  contains  the  estimated  reciprocal  condition  numbers  of  the  selected 
eigenvectors. 

8      STRSEN  Estimating  the  Condition  of  a  Cluster  of  Eigen- 
values 

In  this  section,  we  first  show  the  usage  of  LAPACK  subroutine  STRSEN  for  estimating 
the  reciprocal  condition  number  of  a  specified  multiple  (or  clustered)  eigenvalue  and  its 
corresponding  invariant  subspace,  and  then  outline  the  algorithm. 
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8.1      Usage 
Single  precision 

CALL  STRSEN (SELECT , N . T , LDT . S , SEP , WORK , LDWORK , NWORK , 
$  X.V.ISGN.INFO) 

* 

*  .  .    Scalar  Arguments  . . 

INTEGER  N,  LDT,  LDWORK,  NWORK,  INFO 

REAL  S,  SEP 

* 

*  . .  Array  Arguments  . . 
LOGICAL  SELECTC*) 
INTEGER  ISGN(*) 

REAL  TCLDT,*),  WORK(LDWORK.*) ,  X(*),  V(*) 

* 

*  Arguments 

itn         =  =  =  =  =  =  =  =  = 

* 

*  SELECT  -  LOGICAL  array  if  DIMENSION  (N) 

*  On  entry,  SELECT  specifies  the  1  by  1  or  2  by  2  diagonal 

*  blocks  in  the  eigenvalue  cluster.  For  2  by  2  blocks, 

*  the  first  index  of  SELECT  must  be  set  to  .TRUE,  if  the 

*  block  to  be  collected.  Complex  conjugate  eigenvalues  will 

*  either  both  be  inside  the  cluster,  or  both  outside. 

*  On  exit,  SELECT  may  have  been  altered.  If  the  elements  of 

*  SELECT  corresponding  to  a  2  by  2  block  were  each  initially 

*  set  to  .TRUE. ,  the  program  resets  the  second  one  to  .FALSE. 
* 

*  N      -  INTEGER 

*  On  entry,  N  specifies  the  order  of  the  matrix  T.  N  must  be 

*  at  least  zero. 

*  Not  modified. 
* 

*  T      -  REAL  array  of  DIMENSION(LDT,N) 

*  On  entry,  T  contains  an  upper  quasi-triangular  matrix  in 

*  Schur  canonical  form.  This  means  that  the  diagonal  entries 

*  of  2  by  2  diagonal  blocks  must  be  equal. 

*  Not  modified. 
* 

*  LDT    -  INTEGER 

*  On  entry,  LDT  specifies  the  first  dimension  of  T  as 

*  declared  in  the  calling  (sub)program.  LDT  must  be  at 

*  least  max(l ,  N) . 

*  Not  modified. 
* 

21 


REAL 

On  exit,  S  is  a  lower  bound  on  the  reciprocal  of  the 
condition  number  for  the  selected  cluster  of  eigenvalues. 
S  cannot  underestimate  the  true  reciprocal  condition 
number  by  more  than  a  factor  of  sqrt(N). 

REAL 

On  exit,  SEP  is  the  estimated  reciprocal  of  the  condition 

number  of  the  corresponding  invariant  subspace. 

REAL  array  of  DIMENSION (LDWORK,  N2) 

Workspace,  where  N2  is  the  number  of  eigenvalues  in  the 
cluster,  counting  multiplicities. 


LDWORK  -  INTEGER 

On  entry,  LDWORK  specifies  the  first  dimension  of  WORK  as 
declared  in  the  calling  (sub)program.  LDWORK  must  be  at 
least  rnaxCl ,  N-N2) . 
Not  modified. 

NWORK  -  INTEGER 

On  entry,  NWORK  specifies  the  largest  possible  columns  of 
working  array  U.  NWORK  should  be  larger  than  N2. 
Not  modified. 


*  S 

* 

* 
* 
* 

*  SEP 

*  WORK 

* 
* 
* 

* 
* 

* 

* 
* 
* 
* 

* 


ISGN 


INFO 


REAL 
Workspace 

REAL 
Workspace 

INTEGER 
Workspace 


array  of  DIMENSI0N(N1*N2) 
array  of  DIMENSI0N(N1*N2) 
array  of  DIMENSI0N(N1*N2) 


INTEGER 

On  exit,    INFO   is   set  to 
0  normal  return. 

-K  input  parameter  nxunber  K  is    illegal. 

N+1  there  are  not   enough  columns  for  working 

array  WORK. 


Double  precision.  The  calling  sequence  of  the  double  precision  routine  DTRSEN  is  the 
same  as  that  of  the  corresponding  single  precision  "S"  subroutine  except  that  all  the  real 
variables  are  double  precision. 
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Complex.  The  calling  sequence  of  the  single  precision  complex  routine  is  essentially 
the  same  as  STRSNA,  except  that  the  T,  RE,  WORK,  X,  V  arrays  are  complex.  Integer  array 
ISGN  is  not  needed. 

Double  precision  complex.  The  calling  sequence  of  the  double  precision  comple.x 
routine  ZTRSEN  is  the  same  as  that  of  the  corresponding  single  precision  complex  "C"  sub- 
routine except  that  all  the  real  variables  are  double  precision  and  all  the  complex  variables 
are  double  precision  complex. 

8.2      Example 

The  following  program  segment  illustrates  the  use  of  the  single  precision  subroutine  to 
estimate  the  reciprocal  condition  numbers  of  a  specified  cluster  of  eigenvalues  and  its  cor- 
responding invariant  subspace  for  a  general  matrix.  The  program  first  reduces  a  general 
matrix  to  upper  Hessenberg  form  by  SGEHRD,  and  then  computes  the  Schur  decomposition 
by  the  multishift  QR  algorithm  (SHSEQR).  After  that  the  user  should  input  the  logical  array 
SELECT  to  specify  the  eigenvalues  in  the  cluster.  Then  STRSEN  is  called  to  estimate  the 
reciprocal  condition  numbers. 

PROGRAM        TEST 
INTEGER        ISGN (50) 

REAL  A(S0,50),   WR(50) ,    WI(SO) 

REAL  X(IOO),    V(IOO),   U(50.150),    WORK(50,150) 

LOGICAL        SELECT(50) 

INTEGER        N,    LDA ,    LDWORK,    M,    I,    J,    MAXN.    INFO 
INTEGER        NBLCKl,    NSHIFT,    NBLCK2 
REAL  S,    SEP,    DUMMY 

PARAMETER    (LDA   =   SO,    LDU  =   50,    LDWORK   =   50.    MAXN   =    150) 
* 

*  Input  data: 

*  N:  the  order  of  matrix  A. 

*  A:  N  by  N  array  to  store  the  input  matrix, 
the  blocksize  used  in  Hessenberg  reduction, 
the  number  of  shifts  used  in  the  multishift  QR  algorithm, 
the  blocksize  used  in  bulge  chasing  in  QR  iteration. 


*  NBLCKl 

*  NSHIFT 


*  NBLCK2 
* 

READ(*,*)  N 

DO  10  I  =  l.N 

READ(*,*)  (A(I,J),J  =  l.N) 
10      CONTINUE 

READ(*,*)  NBLCKl 

READ(*,*)  NSHIFT 

READ(*,*)  NBLCK2 
* 

*  Compute  Schur  decomposition. 
CALL  XENVIR( 'BLOCK' ,  NBLCKl) 
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CALL  SGEHRDCH',  N,  A,  LDA,  DUM ,  DUM .  WR.  WORK,  LDWORK , 
$  MAXN,  INFO) 

CALL  XENVIRCSHIFT' ,  NSHIFT) 
CALL  XENVIRC 'BLOCK' ,  NBLCK2) 

CALL  SHSEQRCS',  N,  A,  LDA.  DUM,  DUM,  WR,  WI ,  U.  LDU, 
$  MAXN,  WORK,  LDWORK,  MAXN,  INFO) 

.   DO  20  I  =  1,  N 

WRITEC*,*)  WR(I) ,WI(I) 
20      CONTINUE 

*  Input  the  SELECT  to  specify  the  eigenvalues  to  be  collected 

*  together,  then  the  condition  numbers  of  the  corresponding 

*  invariant  subspace  will  be  estimated. 
* 

DO  30  I  =  l.N 

READ(*,*)  SELECT(I) 
30      CONTINUE 


* 


* 


CALL  STRSEN (SELECT,  N,  A,  LDA,  S,  SEP,  WORK.  LDWORK, 
ISGN.  X.  V.  INFO) 


*      Print  output 

WRITEC*,*)  S.SEP 

STOP 

END 

8.3      Outline  of  the  Algorithm 

STRSEN  routine  estimates  the  reciprocal  condition  numbers  of  specified  multiple  (or  clus- 
tered) eigenvalues  and  their  corresponding  invariant  subspace  for  a  real  Schur  canonical 
matrix  T. 

Logical  array  SELECT  specifies  the  1  by  1  (for  real  eigenvalues)  or  2  by  2  (for  complex 
conjugate  eigenvalues)  diagonal  blocks  that  are  to  be  collected  together  to  form  the  desired 
cluster.  Note  that  for  2  by  2  diagonal  blocks,  the  first  index  of  SELECT  must  be  set  to 
.TRUE,  if  the  block  to  be  collected.  For  real  matrices,  complex  conjugate  eigenpairs  will 
both  be  selected  if  one  is  selected. 

STRSEN  first  calls  subroutine  STREX2  to  collect  the  selected  diagonal  blocks  by  orthogonal 
transformation  to  the  top-left  corner  of  T  such  that 

"^e"  =  ( "^o'  Z 

where  the  selected  blocks  have  been  collected  in  ni  by  nj  matrix  Tn-  Then  STRSEN  computes 
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the  projector  P  on  the  invariant  subspace  associated  with  Tn-  It  is  known  that 

7     -R 


''■00 

where  R  is  the  solution  of  the  Sylvester  equation 

TuR  -  RT22  =  T^2. 

This  is  done  by  subroutine  STRSYL  .  The  program  tests  to  avoid  overflow  if  \\R\\  is  very 
large,  returning  zero  as  the  reciprocal  condition  number. 

The  return  value  S  of  STRSEN  is  the  lower  bound  (1  +  ||i?||/r)~^^^  on  the  reciprocal  of 
||P||2.  It  cannot  underestimate  ||P|I^^  by  more  than  a  factor  of  n^/*. 

Finally,  STRSEN  estimates  the  separation  of  Tu  and  T22-  We  know  that  this  can  be 
estimated  by  the  one-norm  of 

A'~^   =  (/„2   ®  Tn  -  T22  ®  In,  )"\         n  =  Til  +  712 

because  of  the  relationship 

"A'-' 111  <  ,^     ^    ,  =  ||A-^||2  <  v/^M^||A-i||i. 


\/"i"2  sep(Tii,T22) 

This  is  done  by  calling  SLACON  via  a  reverse  communication  interface,  providing  the  solution 
vectors  x  and  y  of  the  equations: 

Kx  =  z,      K   y  =  z 

where  z  is  determined  by  SLACON.  This  means  we  must  solve  the  Sylvester  equations: 

TiiA'  -  A'r22  =  Z     T^^Y  -  YTJ2  =  Z 

This  is  again  is  done  by  the  subroutine  STRSYL. 

The  return  value  SEP  of  STRSEN  is  the  estimated  (upper  bound)  of  sep(Tii,r22). 

Swapping  adjacent  diagonal  blocks  on  the  diagonal  of  the  input  matrix  costs  0{n),  so 
swapping  Ui  selected  eigenvalues  to  the  top  left  corner  costs  at  most  O(nin^)  (and  as  little 
as  nothing  if  they  are  already  at  the  top  left  corner).  Once  the  desired  eigenvalues  are  at  the 
top  left,  solving  either  above  Sylvester  equation  costs  0(nj7i2  + 711x12)  operations.  Therefore 
the  condition  number  estimation  of  a  cluster  of  eigenvalues  and  their  corresponding  invariant 
subspace  costs  at  most  0(71'')  operations,  or  as  few  as  O(n^)  if  7ii  C  7i2. 
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A      Solution  of  the  Sylvester  Equation 

Considerable  attention  in  tlie  literature  has  been  paid  to  solving  the  Sylvester  equation. 
Among  proposed  solutions,  Bartels  and  Stewart's  method[2],  and  Golub,  Nash  and  Yaw 
Loan's  method[r2]  are  direct  matrix  factorization  methods.  In  this  appendix,  we  discuss 
the  method  originally  presented  by  Bartels  and  Stewart,  and  the  associated  routine  STRSYL. 
The  Sylvester  matrix  equation  is  of  the  form 

op(A)X  -  Xop(B)  =  sC  (A.IO) 

where  A  ,  B  and  C  are  real  m  x  m,  7i  x  n  and  m  x  n  matrices  respectively.  op(A)  =  A 
or  A'^  is  a  transpose  option.  5  is  a  scaling  factor  (<  1)  which  is  so  chosen  so  that  A'  can 
be  computed  without  overflow.  We  will  also  suppose  that  A  and  B  are  in  upper  quasi- 
triangular  form,  otherwise  we  should  compute  the  Schur  decomposition  of  .4  and  B  (by 
SGEHRD  and  SHSEQR), 

U^AU  =  R,      V'^BV  =  S  (A. 11) 

where  R  and  S  are  upper  quasi-triangular,  and  U  and  V  are  orthogonal.  The  reductions 
(A.ll)  lead  to  the  system  (A.IO). 

It  is  well  known  that  (A.IO)  has  a  unique  solution  if  and  only  if  there  are  no  common 
eigenvalues  of  A  and  B. 

Let  p  be  the  number  of  1  by  1  and  2  by  2  blocks  along  the  diagonal  of  A,  and  let  q  be 
the  number  of  the  1  by  1  and  2  by  2  blocks  along  the  diagonal  of  B.  Partition  the  A,  B,  X 
and  C  conformally.  If  op(yl)  =  A  and  op{B)  =  B  then  the  ijth  block  A',j  satisfies 

^., A',,  -  X,,B,,  =  51  (C„  -  C'ij)  ( A.12) 

where 

p  i-i 

C;^=    XI    A,,X,,-Y,XaBi, 

k=t+l  1=1 

Note  that  since  ^„  and  Bjj  are  each  1  by  1  or  2  by  2,  the  solution  of  (3)  can  be  obtained  by 
solving  a  linear  system  of  order  at  most  four.  That  can  be  solved  easily.  Once  calculated,  the 
solution  Xij  can  be  stored  in  the  locations  occupied  by  C,j,  which  is  no  longer  needed.  The 
solution  matrix  A'  can  be  successively  solved  column  by  column  starting  from  bottom-left 
corner  of  A',  i.e.,  in  order  Xpi,  A'p_i,i,  . . . ,  A'n,  A'p2,  . . . ,  A'l,. 

For  solving  equation  (A.12),  we  note  that  if  >!„  and/or  Bjj  are  balanced  2  by  2  blocks, 
i.e.,  they  are  of  the  form 

'  Q     /3 

7      Q 

then  a±  n  are  the  eigenvalues  where  0f  =  -n^.  Equation  (A.12)  can  then  be  expressed  as 
a  2  by  2  linear  system 

{H  -  w)Y  =  SiF 

Here  //  is  an  na  by  na  real  matrix  {na  =  1,2),  w  is  real  or  complex,  Y  and  F  are  na  by  1 
matrices  which  are  real  if  w  is  real  and  complex  if  w  is  complex,  and  5i  is  a  local  scaling 
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factor  (<  1)  which  is  so  chosen  that  Y  can  be  computed  without  overflow  (see  SLALN2  and 
SLASY2).  In  particular,  if  A,-,  and  Bjj  have  the  same  eigenvalues,  Si  is  set  to  0. 

Similarly,  if  op(.4)  =  .4'^  and  op(5)  =  B,  the  ijtlx  block  A',j  of  the  solution  A'  can  be 
successively  solved  column  by  column  starting  from  top-left  corner  of  X ,  i.e.,  in  order  A'n- 
A  21 ,  . . . ,  A  pi ,  A 1 2 ,  . . . ,  A  p^ . 

If  op(/l)  =  A  and  op(i?)  =  B^ ,  the  ijth  block  A',j  of  the  solution  A'  can  be  successively 
solved  column  by  column  starting  from  bottom-right  corner  of  A',  i.e.,  in  order  A'p,,  A'p_i,g, 

.  .  .  ,  Aiq,  Ap_^_i,   .  .  .  ,  All. 

If  op(/l)  =  A^  and  op{B)  —  B^ ,  the  ijtli  block  A',j  of  the  solution  A'  can  be  successively 
solved  column  by  column  starting  from  top-right  corner  of  A',  i.e.,  in  order  A'l,,  A'2,,  . . ., 

^pqi  -^1,9-11  •  •  •  7  -^pl- 

Using  the  above  different  substitution  orderings  enables  one  to  work  on  the  input  ma- 
trices directly  rather  than  to  transpose  the  input  matrices. 

The  overall  number  of  flops  for  the  above  substitution  solution  is 

where  we  have  assumed  that  A  and  B  are  already  in  Schur  form. 

The  program  may  be  used  to  iteratively  refine  of  the  computed  solution  A'l  of  (A. 10):  let 
the  residual  matrix  R\  —  C  —  AXi  +  XiB  be  computed  in  double  precision  and  rerounded 
to  single  precision.  Use  the  same  program  to  solve  the  system  AA'2  —  X2B  —  Ri.  Then 
Xi  +  X2  will  in  general  be  a  more  accurate  approximation.  This  process  may  be  iterated. 
This  iteration  is  analogous  to  the  iterative  refinement  of  approximate  solutions  of  linear 
system  as  described  by  Wilkinson[27,  p.  255].  (This  is  not  done  in  STRSEN  and  STRSNA.) 

B      Swapping  Diagonal  Blocks 

The  crux  of  swapping  a  selected  block  of  a  real  Schur  form  to  a  specified  position  along  the 
diagonal  (subroutine  STREXC),  or  collecting  selected  blocks  together  (subroutine  STREX2) 
is  the  swapping  of  adjacent  blocks  by  an  orthogonal  similarity  transformation  (subroutine 
SLAEXC).  Stewart[23]  developed  an  adjacent  block  swapping  algorithm  using  one  or  two  QR 
steps  with  a  pre-determined  shift  to  force  the  ordering  of  the  eigenvalues  of  the  new  blocks. 
More  recently,  Ng  and  Parlett[17]  present  a  more  straightforward  algorithm  for  the  same 
task.  The  presentation  in  this  appendix  is  based  on  Ng  and  Parlett's  work.  We  discuss  in 
more  detail  the  treatment  of  pathological  cases. 
Consider  a  submatrix  of  the  form 

T\\     T\2 
0      T22 

where  T\\  is  a  p  by  p  matrix,  and  r22  is  a  g  by  q  matrix,  p,  9  =  1  or  2,  and  assume  that 
Til  and  T22  have  no  eigenvalue  in  common.  Moreover,  we  assume  that  if  either  is  a  2  by  2 
matrix,  it  has  been  standardized  (i.e.,  it  has  identical  diagonal  entries.)  Now,  we  want  to 
find  an  orthogonal  matrix  Q  which  swaps  T\\  and  T22,  i.e.. 


where  T„  is  similar  to  T,;,  i  =  1,2. 

Since  Tn  and  T22  have  no  eigenvalue  in  common,  it  follows  that  there  exists  a  unique 
;;  X  q  matrix  A'  such  that 

TiiA  -  AT22  =  Ti2- 


Hence 


^11     T,2\     ^     /   Ip     -X   \       Tu       0     W  7p     A 
0      T22  I  V    0       I,     I  [    0      T22  )  [   0     I, 

-A     Ip\  (  T22      0     W   0     /, 


/,       0    ;  \     0      Tn    j  \  Ip    A 
We  see  that  it  is  easy  to  find  an  orthogonal  [p  +  q)  x  [p  +  q)  matrix  Q  such  that 

with  some  invertible  qxq  M2,  e.g.,  using  Householder  matrices  to  do  the  QR  decomposition. 
Let  Q  premultiply  and  postmultiply  the  original  matrix,  yielding 

A'     Ip\  (  T22      0    W   0     /, 


71/2    w  \  (  T22    0   W  AI2    w 

0      A'h  )  [    0      Tn  )[    Q      Ml 

Mn     W   \  (  T22      0    \  (  M^'     -M.-^WM-^ 
0      M,  )  [    0      Tu  )[      0  M-i 

M2T22M^^  T[. 

KhTnMi' 

Til  '^nd  T22  have  been  swapped. 

The  above  considerations  are  summed  up  in  the  following  steps. 

1.  Solve  TiiA'  —  XT22  =  •sTi2.  5  is  a  scale  factor  introduced  to  avoid  overflow. 

2.  Check  if  the  magnitude  of  ||A'||  exceeds  the  square  root  of  the  overflow  threshold.  \\\ 
this  case  T\\  and  T22  are  too  close  to  swap,  so  we  exit. 

3.  Use  a  Householder  matrix  Q  to  do  the  QR  decomposition  of  (A'    I)^  and  update  T 
by  Q:  QTQ^, 

4.  Accumulate  the  orthogonal  transformations  if  desired. 

5.  To  preserve  the  standard  Schur  form,  make  the  diagonal  elements  equal  in  each  2  by 
2  block  using  orthogonal  transformations. 

6.  Accumulate  the  orthogonal  transformations  if  desired. 
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Several  comments  should  be  made.  First,  the  solution  of  the  matrix  equation  T^X  - 
XT22  =  sTi2  has  been  discussed  in  detail  in  Appendix  A,  the  routine  SLASY2.  Second. 
there  is  no  danger  in  working  with  A'  of  large  norm  provided  that  ||A'||"  does  not  overflow. 
Moreover  if  HA'!!'  does  overflow  then  the  blocks  should  not  be  swapped  because  a  tinv 
perturbation  will  cause  Tn  and  T22  to  have  at  least  one  common  eigenvalue[9].  Hence  in 
step  2,  we  check  the  norm  of  A',  and  if  A'  satisfies 

5-max(||rn|M|T22||) 

l|A'||  +  5  ^'' 

where  e  is  the  machine  precision,  then  the  two  blocks  are  regarded  as  too  close  to  swap. 
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C      List  of  LAPACK  Routines  for  the  Nonsymmetric  Eigen- 
problem 

LAPACK  main  routines  for  the  nonsymmetric  elgenproblem: 

SGEBAL  Balance  an  input  general  matrix  and  isolate  eigenvalues  whenever  possible. 

SGEBAK  Form  the  eigenvectors  for  a  general  matrix  by  back  transforming  those  of  the  cor- 
responding balanced  matrix  determined  by  SGEBAL. 

SGEHRD  Reduce  a  general  matrix  to  an  upper  Hessenberg  matrix. 

SHSEQR  Compute  the  eigenvalues  of  an  upper  Hessenberg  matrix  by  the  multishift  QR 
algorithm,  and  return  the  Schur  form,  accumulating  the  orthogonal  matrix  if  desired. 

STREVC  Compute  selected  right  and/or  left  eigenvectors  of  a  Schur  matrix. 

SHSEIN  Compute  selected  right  and/or  left  eigenvectors  of  a  Hessenberg  matrix  by  inverse 
iteration. 

S0RGC3  Overwrite  a  matrix  containing  Householder  vectors  stored  in  the  strictly  lower  part 
by  the  orthogonal  matrix  they  represent. 

STRSNA  Estimate  selected  reciprocal  condition  numbers  of  individual  eigenpairs  of  Schur 
matrix. 

STRSEN  Estimate  selected  reciprocal  condition  numbers  of  a  multiple  (or  cluster  of)  eigen- 
values and  their  corresponding  invariant  subspace  of  a  Schur  matrix. 

STRSYL  Solve  the  Sylvester  equation  with  coefficient  matrices  in  Schur  form. 

STREXC  Swap  a  selected  diagonal  1  by  1  or  2  by  2  block  of  a  Schur  matrix  to  a  specified 
position. 

STREX2  Collect  several  selected  diagonal  1  by  1  or  2  by  2  blocks  of  a  Schur  matrix  to  the 
top-left  or  bottom  right  corner. 

LAPACK  auxiliary  routines  for  the  nonsymmetric  eigenproblem: 

SLAHRD  Chase  a  k  by  k  bulge  of  an  upper  Hessenberg  matrix  one  block  down  from  a  specified 
column  number. 

SLAHQR  BLAS   1  based  QR  routine  to  compute  the  eigenvalues  of  an  upper  Hessenberg 
matrix,  and  return  the  Schur  form,  accumulating  the  orthogonal  matrix  if  desired. 

SLATRS  Mixed  subroutine  of  BLAS  1  and  BLAS  2  to  solve  triangular  equations  while  avoid- 
ing overflow. 

SLAQTR  Solving  real  or  complex  quasi-triangular  systems  where  the  real  part  is  quasi- 
triangular,  and  the  imaginary  part  is  of  a  special  form. 
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SLALN2  Solve  a  2  by  2  linear  equation. 

SLAE2  Compute  the  eigenvalues  of  a  2  by  2  nonsymmetric  matrix. 

SLAEXC  Swap  adjacent  diagonal  1  by  1  or  2  by  2  blocks  of  a  Schur  matrix. 

SLASY2  Solve  the  Sylvester  equation  with  coefficient  matrices  up  to  2  by  2. 

SLAEQU  Equalize  the  diagonal  elements  of  a  2  by  2  block  with  an  orthogonal  similarity. 

Other  LAPACK  routines,  auxiliary  routines,  functions  called  by  eigensysteni 
subroutines  (except  Level  1,  2  or  3  BLAS  routines). 

XERBLA,    LSAME,    RIMACH,    ENVIR 

SLACON,    SLACPY,    SLAZRO .    SLARFG.    SLARF ,    SLANHS ,    SLAPY2 ,    SLAPY3 
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